Scaling and self-similarity in an unforced flow of inviscid fluid trapped inside a viscous 

fluid in a Hele-Shaw cell 
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We investigate quasi-two-dimensional relaxation, by surface tension, of a long straight stripe of 
inviscid fluid trapped inside a viscous fluid in a Hele-Shaw cell. Combining analytical and numerical 
solutions, we describe the emergence of a self-similar dumbbell shape and find non-trivial dynamic 
exponents that characterize scaling behavior of the dumbbell dimensions. 
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Introduction. Consider a bubble of low-viscosity fluid 
(say, water) trapped inside a high-viscosity fluid (say, 
oil) in a quasi-two-dimensional Hele-Shaw cell. What 
will happen to the shape of the bubble, if the (hori- 
zontal) plates are perfectly smooth, and the fluids are 
immiscible? The answer depends on the initial bubble 
shape. A perfectly circular bubble (or an infinite straight 
stripe) will not change, while a bubble of any other shape 
will undergo surface-tension-driven relaxation until it ei- 
ther becomes a perfect circle, or breaks into two or more 
bubbles, which then become perfect circles. The bubble 
shape relaxation is non-local, as it is mediated by a vis- 
cous flow in the outer fluid. The resulting free boundary 
problem is hard for analysis. This is especially true when 
the bubble has a complex (even fractal) shape, like that 
observed, in radial geometry, in a strongly forced Hele- 
Shaw flow, when the viscous fluid was initially displaced 
by the inviscid fluid Q. The shape complexity results 
from the viscous fingering instability 0, 01 ■ The forced 
Hele-Shaw flow is a celebrated problem in fluid dynamics 
and nonlinear dynamics 0, 0, IE • The role of small 
surface tension there is to introduce a (nontrivial) reg- 
ularization on small scales. This Letter deals with an 
unforced Hele-Shaw (UHS) problem, where surface ten- 
sion is the only driving mechanism. We address the UHS 
problem in the case when the inviscid fluid is initially in 
the form of a long stripe. We show that this special initial 
condition provides a useful characterization of the UHS 
model, as the evolving stripe, which develops a dumbbell 
shape, exhibits self-similarity with non-trivial dynamic 
exponents. 

UHS problem. Let the inner fluid have negligible vis- 
cosity, so that the pressure inside the bubble is ho- 
mogeneous. The velocity of the viscous outer fluid is 
v(r, t) = — (b 2 /\2y) Vp (r, i), where p is the pressure, 
is_ the_dynamic viscosity, and b is the plate spacing 
Therefore, the interface speed is 



-(&7l2/x)V„p, 



(1) 



of the outer fluid, the pressure is a harmonic function: 

V 2 p = . (2) 
The Gibbs-Thomson relation at the interface yields 



p | 7 = (it/4) o-K, . 



(3) 



where a is surface tension, and K. is the local curvature of 
the interface, positive when the inviscid region is convex 
outwards. As both the supply of the inner fluid, and 
evacuation of the outer fluid are blocked, we demand 



V„p | r = 



(4) 



where index n denotes the components of the vectors nor- 
mal to the interface, and V„p is evaluated at the respec- 
tive points of the interface 7. In view of incompressibility 



at the external boundary of the system L. Equations 
Q-Q define the exterior UHS problem (see Ref. || 
for a more detailed discussion). A related, but different 
interior problem has been also considered, mainly in the 
context of singularity formation (pinch-offs) in bubbles 
of viscous fluid |9( . The UHS model has two important 
properties: (i) the bubble area remains constant, (ii) the 
length of the interface between the two fluids is a non- 
increasing function of time |ic|. 

The UHS problem is not integrable. Moreover, we are 
unaware of any analytical solutions to this problem, ex- 
cept for a linear analysis of a slightly deformed flat or 
circular interface Owing to its two-dimensionality, 

the problem can be reformulated as a nonlocal nonlinear 
partial differential equation for a conformal map which is 
analytic in the exterior of the unit circle [l2] | . This equa- 
tion, however, is hard for analysis. We consider here a 
simple but non-trivial case that can be analyzed directly 
in the physical plane: the dynamics of a half-infinite (or, 
physically, very long) stripe. 

Stripe dynamics: theoretical predictions. Let at t = 
the bubble have the form of a half-infinite straight 
stripe of width A, located along the rc-axis as shown 
in Fig. 2] The external boundary of the system T is 
at infinity, where the pressure is bounded. We will 
measure the distance in units of A, the time in units 
of t — 48/iA 3 /(7rcr& 2 ), and the pressure in units of 
Po = 7rer/(4A). In the rescaled variables Eqs. Q and 
© become v n = —V n p and p | 7 = JC, so the rescaled 
problem is parameter-free. 
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FIG. 1: Setting for the dumbbell dynamics. 



We are interested in the late-time behavior: i > 1. 
Because of the Gibbs-Thomson effect, the pressure gra- 
dient is largest near the tip, so the tip retreats along the 
x-axis. As the bubble area must be conserved, the re- 
treating stripe acquires a dumbbell shape, and the lobe 
of the dumbbell expands with time, see Fig. |21 Surpris- 
ingly, the main contribution to the dumbbell area comes, 
at t ^> 1, from the dumbbell neck, and not from the lobe. 

Going over to a quantitative analysis, we assume (and 
later verify numerically) that the lobe can be character- 
ized by a single time-dependent length scale R(t). An- 
other length scale is L(t): the retreat distance of the 
dumbbell. Our main objective is to find the exponents 
of the power laws for R(t) and L(t). Our analysis will 
not give the numerical coefficients of these power laws 
(which, in the rescaled units, are of order unity); these 
will be found numerically. Introduce polar coordinates 
r and (f>, see Fig. ^ The dumbbell neck, r ^> R(t), is 
almost flat, so p must vanish at — > and (f> — > 2ir. On 
the other hand, p — JC ~ 1/R(t) at the lobe interface (for 
definiteness, at <f> = ±7t/2). Therefore, the leading term 
in the multipole expansion 13] is 



p{r^,t) = C[R{t)ry 1/2 sin(0/2), 



(5) 



where C — 0(1). Having demanded the Gibbs-Thomson 
condition here, we somewhat stretched the validity of Eq. 
©, but this can only affect the value of constant C . The 
dashed lines in Fig. ^ show the field lines of Vp. 

Equation JSJ yields the normal component of the inter- 
face speed v n = —V n p in the neck region. For the upper 
interface of the neck 

U " = '^ ( ^° ) = ' 2^4 )r 3/ 2 - (6) 

Now we return to the Cartesian coordinates. Let h(x±, t) 
be the local height of the dumbbell, while xi = x — L(t) 
be the horizontal coordinate in the moving frame with 
the origin at the tip. In the neck region, x\ 3> R(t), the 
quantity dh(xi, t)/dt is given by Eq. ©, so we obtain 



h(xut) -- = 



Colt' 



2R 1 / 2 (t' 



3/2 



(7) 



This equation yields h{x\,t) in two different limits. At 
very large distances from the lobe, x 3> L(t) (region I) 



h{x] 



' 2 



C 



2x 



3/2 



dt' 



3/2 



R}/ 2 {t) 



(8) 



where the last estimate assumes that R(t) is a power of 
t. Another limit corresponds to intermediate distances: 
R(t) <C x\ <C L(t) (region II). Here, at fixed x, the main 
contribution to the integral in Eq. JJJ comes from times 
close to t, so that X\{t)/L{t) <C t — t' <C t . Indeed, one 
can expand x\(t') = x\{t) + L(i)(t — <') + ... and, in 
the leading order, ignore higher order terms. The effec- 
tive time interval for the integration is (t — St' , t), where 
St' - Xx(t)/L(t). Furthermore, R 1 / 2 ^') can be evaluated 
at t' — t, as its variation on the time interval (t — 5t',t) 
is negligible. Then, extending the lower limit of the in- 
tegral to — oo and calculating the remaining elementary 
integral, we obtain 



C 



R l / 2 {t)L(t)x\ ,2 {t) 



(9) 



Now we can estimate the contributions of regions I and 
II to the dumbbell area gain A in the neck region. We 
integrate Eq. ((HJ over x\ from, say, 2L(t) to infinity, and 
Eq. © from R(t) to 2L(t). The results are: 



and 



A n (t) 



LV2(t) fll/2(t) 



L^jt) 
L(t) Ry 2 (t) 



in region I , (10) 



in region II . 



(11) 



Once L(t) is a power law, Aj and Aji are comparable. 
Notice that in region I (respectively, II) the main contri- 
bution comes from the lower (respectively, upper) limit 
of integration. As we verify a posteriori, the contribution 
to the dumbbell area of the lobe itself, An ~ R 2 (t), is 
negligible compared to Ai and An as long as t ^> 1. 

Now we can find the dynamic exponents of L{t) and 
R(t). First, we employ the area conservation of the 
dumbbell. The area loss L(t) x 1 of the retreating dumb- 
bell must be equal to the area gain in the neck, so up to 
numerical coefficients of order unity 



L(t) - Aj(t) - A H (t) 



LV2(t) flV2(t) 



(12) 



Second, there is a simple kinematic relation between L(t) 
and the characteristic speed of the lobe motion V/. Us- 
ing Eq. ||SJ|, we obtain V/ ~ —dp/dr (r ~ R(t), 4> — 7r) ~ 
R~ 2 (t), and demand L(t) ~ R~ 2 (t). Combined with 
Eq. ijnj, this yields 



L{t) - t 3/5 and R(t) - t L/b at t > 1 



1/5 



(13) 
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Now we can return to Eqs. (|5|)- (|llfl and find the explicit 
time-dependences. For example, the far-neck asymptote 
in Eq. © becomes h(x,t) - 1/2 ~ t 9 / 10 x- 3 / 2 . We can 
also verify that, at t ^> 1, the lobe area Ar ~ R 2 {t) ~ 
t 2/5 is indeed much less than Aj(t) ~ A H (t) ~ i 3/5 . 

That the lobe is characterized by a single dynamic 
length scale R(t) ~ i 1 / 5 implies a similarity Ansatz for 
the lobe shape in the moving frame: 

h a (xi,t)=t 1 / 5 $(x 1 /t 1 / 5 ) ati>l. (14) 

Numerical method. To test our predictions, we per- 
formed simulations of the dynamics of long stripes with 
dimensions X x 1, where X ^> 1. The ultimate shape 
of such a stripe is a perfect circle. Therefore, the scaling 
behavior, predicted by our theory of a one-sided dumb- 
bell, appears as an intermediate asymptote, as we require 
R(t) ^> 1 but L(t) <C X. In view of the predicted scalings 
with time, we must demand 1 <C t <C X 5 / 3 . 

Our numerical algorithm 0] is based on a representa- 
tion of the harmonic field in terms of a line integral over 
the bounding contour; it involves tracking of the contour 
nodes. We employed a variant of the boundary integral 
method, suggested in Ref. ^^|. The algorithm includes 
solving an integral equation for an effective density of 
dipole moments (DMD) and evaluating another integral, 
which yields a harmonic conjugated function (HCF). The 
normal velocity of the interface is given by the derivative 
of the HCF along the contour. The very large aspect 
ratio of the dumbbell demands a different discretization 
compared to Ref. 0. Indeed, the typical scale of vari- 
ation of the kernel of the integral equation 0] over the 
almost flat neck of the dumbbell is close to 1: the initial 
stripe thickness. On the other hand, the DMD changes 
much slower there. This enabled us to considerably re- 
duce the number of grid nodes in the neck region. Wc 
used a piecewise constant function to approximate the 
DMD, and a piecewise linear function to approximate 
the contour. Therefore, each of the integrals was approx- 
imated as a discrete sum of the DMD values multiplied by 
an integral of the kernel between two neighboring nodes. 
The latter integrals can be calculated analytically. The 
HCF is evaluated at middle points between the nodes, 
while the normal velocity at each node is evaluated using 
the values of the HCF at neighboring middle points. 

We used an explicit finite difference method to track 
the contour. The number of grid points, needed for an 
accurate solution and contour tracking, decreases with 
time together with the perimeter of the dumbbell. An 
obvious modification of the algorithm of Ref. 0] ex- 
ploited the 4-fold symmetry of the dumbbell. The area 
conservation of the dumbbell was used for accuracy con- 
trol. The time step chosen was 5 x 10 -3 min \R l /v % n \, 
where Ri and v l n are the local curvature radius and nor- 
mal velocity, respectively, in the node i of the contour. 
This choice resulted in good area conservation: in the 
simulation described below less than 0.5% of the area 
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FIG. 2: Snapshots of a part of the simulated system at t = 
0, 1000 and 3010. Notice the large difference between the 
horizontal and vertical scales. The inset shows, to scale, the 
lobe region at t = 3010. 
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FIG. 3: Figure a shows, in a log-log scale, the retreat dis- 
tance L versus time and its power-law fit 2.75 1 ' 60 . Figure b 
shows, in a log-log scale, the maximum dumbbell height, h max 
(the empty circles), and the position of the maximum, x™ ax 
(the filled circles), versus time, as well as their power-law fits 
0.66 1 ' 21 and 0.94 1 - 20 , respectively. 

was lost by the time t = 7000. As the dumbbell contour 
becomes smoother, the time step greatly increases. 

Numerical results. Here we report a simulation with 
X = 2000. Figure |21 shows snapshots of a part of the sys- 
tem at times t = 0, 1000 and 3010. The stripe develops a 
dumbbell shape (though some may prefer a comparison 
with daisy petal). The lobe grows with time, the neck 
widens. Shown in Fig. |2t is the retreat distance L(t) 
versus time. A power law fit yields exponent 0.60 which 
coincides with the theoretical value 3/5. FigureEb shows 
the maximum dumbbell height, h max , and the position 
of the maximum, x™ ax , versus time. Each of these two 
quantities exhibits a power law; the fitted exponents are 
0.21 (for h' max ) and 0.20 (for x? ax ) 7 in agreement with 
the theoretical value 1/5. At long times, when the as- 
pect ratio of the dumbbell is already not large enough, 
the straight line in Fig. |2K slightly curves down, while 
those in Fig. |3t> curve up. We verified that for a shorter 
stripe, X — 1000, deviations from the same straight lines 
occur earlier, as expected. The time interval of the three 
fits, 20 < t < 1000, corresponds to the common parts of 
the dependences for the two values of X. 

Figure 01 depicts the (rescaled) dumbbell shape in the 
moving frame at three different times. The collapse of 
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FIG. 4: Self-similarity of the lobe. Shown is the shape func- 
tion h(xi,t), rescaled to the maximum dumbbell elevation, 
versus the coordinate xi, rescaled to the abscissa of the max- 
imum, at times 160.3 (the filled circles), 1000 (the squares), 
and 3010 (the empty circles). 
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FIG. 5: The dumbbell neck shape and dynamics. Figure a: 
the neck shape at t = 13.9, found numerically (the dotted line 
with circles), Eq. I15H (the solid line), and the first term of 
Eq. I|15|l (the dashed line). Figure b and its inset show, in a 
log-log and linear scales, respectively, the pre-factors a and b 
versus time (symbols). Also shown are a power-law fit with 
exponent 0.92 (figure b) and a linear fit (inset). 



three different curves into a single one supports the sim- 
ilarity Ansatz l|14|). Notice that the dumbbell shape 
h(xi,t) in region II [R(t) <C xi <C L(t)] belongs to the 
similarity region. Indeed, assuming that - £~ 1/2 at 
£ 3> 1, we see that Eq. ifHf) yields Eq. © (where one 
should substitute L ~ t 3 / 5 and R ~ t 1 / 5 , and neglect 1/2 
in the left hand side). 

The self-similarity breaks down at a distance Xi ~ 
L(t) ~ t z / b from the tip. Beyond this distance, Eq. © 
predicts a power-law neck shape. Figure [SJi shows the 
shape of the dumbbell neck at time t — 13.9, computed 
numerically Also shown are the quantity 



-3/2 



(X-2L- x x 



-3/2 



b [xx 2 + (X - 2L - Xl y 2 } , (15) 



and its first term, proportional to a. Equation 115fl differs 
from Eq. JSJ in that (i) it accounts for contributions 
from two dumbbell lobes, and (ii) it accounts for the sub- 
leading term in the multipole expansion of the harmonic 
function p. Note that a and b are the only adjustable 
parameters in Eq. i|15|) . The resulting profile is almost 
indistinguishable from the numerical profile. The first 



term of Eq. (|15(l already gives fairly good agreement. 
The excellent agreement holds, on a shrinking interval 
of xi, until t = 7020. Figure |3Jd and its inset show a 
and b versus time, respectively. A power-law fit of a(t) 
yields exponent 0.92, close to our prediction 9/10. The 
pre-factor b(t) behaves linearly with time. How does it 
compare with the theory? As p ~ 1/R(t) at the lobe 
interface, the sub-leading term p ~ r -1 sin^> does not 
include R(t). Then, repeating the procedure which led 
us to Eq. JHJ, we do obtain b(t) ~ t. 

Summary. We studied the UHS flow in the case when 
the inviscid fluid is initially in the form of a long stripe. 
We found that the resulting dumbbell dynamics exhibit 
self-similarity with nontrivial exponents. The solution 
we obtained is the first analytical solution for an UHS 
flow that goes beyond a linear theory. Similarly to other 
curve-shortening area-preserving relaxation models |lfij . 
the stripe relaxation provides a useful characterization 
of this non-integrable flow. Its experimental realization 
should not be difficult. 
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